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ABSTRACT 

We present an extension to the short-characteristic ray-tracing and non-equilibrium photon- 
ionization code C 2 -Ray. The new version includes the effects of helium and improved multi- 
frequency heating. The motivation for this work is to be able to deal with harder ionizing 
spectra, such as for example from quasar-like sources during cosmic reionization. We review 
the basic algorithmic ingredients of C 2 -Ray before describing the changes implemented, 
which include a treatment of the full on the spot (OTS) approximation, secondary ioniza- 
tion, and multi-frequency photo-ionization and heating. We performed a series of tests against 
equilibrium solutions from CLOUDY as well as c omparison s to the hydrogen only solutions 
by C 2 -Ray in the extensive code comparison in llliev et ail (120061) , We show that the full, 
coupled OTS approximation is more accurate than the simplified, uncoupled one. We find 
that also with helium and a multi-frequency set up, long timesteps (up to ~ 10% of the re- 
combination time) still give accurate results for the ionization fractions. On the other hand, 
accurate results for the temperature set strong constrains on the timestep. The details of these 
constraints depend however on the optical depth of the cells. We use the new version of the 
code to confirm that the assumption made in many reionization simulations, namely that he- 
lium is singly ionized everywhere were hydrogen is, is indeed valid when the sources have 
stellar-like spectra. 
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1 INTRODUCTION 

Photo-ionization is one of the major radiative feedback processes 
in astrophysics. The extreme ultraviolet (EUV) photons produced 
by massive stars in star formation regions are capable of heating 
the gas to temperatures around 10 4 K and the heated ions pro- 
duce copious amounts of collisionally excited line radiation, lead- 
ing to the well-known and sometimes spectac ular images o f H II 
regions (such as for example the Orion nebula. lO'DellluOOlh . Ac- 
creting black holes and neutron stars also produce ionizing radi- 
ation and, depending on their mass, can ionize smaller or larger 
regions around themselves. On galactic scales, the emission from 
supermassive central black holes are observed to produce ioniza- 
tion cones stretching into the galaxy's immediate environment (see 
e.g. |Poggefl 989). Even at the largest scales photo-ionization is im- 
portant. Some time before redshift 6, ionizing radiation escaped 
from the first generations of galaxies and percolated through the 
intergalactic medium (IGM), changing it from cold and neutral to 
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warm and ionized. This process is known as reionization and was 
the Uni verse's last glo bal phase transition (see for example the re- 
view in Barkana 2006). 

Traditionally photo-ionization codes concentrat ed on solving 
equili brium situations , as for example CLOUDY jFerland et alj 
Il998l) . MAPPINGS ( [Sutherland & Dopital 1 19931) and the three- 
dimensional code MOCASSIN tercolano et alj 120081) do. The 
main aim of these codes is to accurately calculate line strengths for 
comparison to spectroscopic observations. However, since the in- 
crease in pressure can drive powerful flows in the gas, there is also 
a need to couple photo-ionization calculations to hydrodynamics. 
This necessitates a simpler version of the radiative transfer, since 
it has to be calculated in step with the hydrodynamics and for dy- 
namic calculations the individual line strengths are less interesting. 
The history of these type s of calculations goes back quite far, see 
for example lYorkel i ll 9861) for an overview of the work done before 
1986. 

For some of the applications mentioned above, a lower dimen- 
sional approach (one or two dimensional) is sufficient. However, 
other applications, such as cosmic reionization, require the trans- 
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port to be performed in the full three dimensions. Because of the 
higher dimensionality here as well a simpler version of the radiation 
and photo-ionization physics is typically implemented, although 
the level of sophisticatio n varies between methods ( see e.g. the first 
generation methods of | Razoumov & Scotj 19991 ; ISokasian et alj 
l200lLlNakamoto et alj|200ll ; lciardi et alj|200lh . In some cases the 
coupling t o the dynamics was also done for cosmological applica - 
tions (e.g.. lRicotti et alj2002l : lTrac et al.l200l ; lwise & Abell2008l) . 

One of the simplifications which is often employed is to 
consider hydrogen as the only element being photo-ionized. 
Since close to 10% of the gas is helium, this approxima- 
tion is crude. However, as shown for example in figure 2.4 of 
lOsterbrock & Ferlandl d2006t) . for typical O-star spectra the ioniza- 
tion of helium follows largely that of hydrogen. So, if one is only 
interested in the structure of the ionized regions, assuming that he- 
lium follows hydrogen is not a bad approximation. 

However, as one moves to harder spectra, this assumption be- 
comes less and less valid. Not only does one need to take into ac- 
count that helium can become doubly ionized, also the fact that 
the cross section of the two types of helium start contributing sub- 
stantially to the opacity of the gas becomes an issue. This prob- 
lem becomes especially important when the ionizing spectrum is 
powerlaw-like, as one expects from hot accretion disks around 
black holes. Specifically for the case of cosmic reionization, where 
there is a possible contribution from powerlaw-like sources such 
as quasars and mini-quasars, any proper study of their contribution 
should consider both hydrogen and helium. 

This has motivated us to extend the capabilities of the code 
C 2 -RAY to include the effects of helium. C 2 -RAY is a photon- 
conserving radiative transfer code that u ses short characteristic 
ray tracing and is described in detail in iMellema et ail (2006b) 
(hereafter M+06 ) . It has been used extensively for reioniza tion 
simulations (e.g. IMellema et alj|2006cl ; llliev et alj|200l 120081 to 
name a few). It was also combined w ith a hydrodynamics code 
(CAPREOLE-C 2 , IMellema et alj|2006al). and an ideal magnetohy- 
drody namics code (Phab-C 2 , Ide Colle & Ragal2006l : lArthur et alj 
1201 1), to investigate galactic H II regions. Furthermore, it was 
tested aga i nst ot her non-equilibrium radiative transfer codes in 
llliev et alj (2006) (hereafter 1+06) and in conjunction with a grid- 
based hydrodynamic code (for details, see lMellema et a l. 2006a) in 
a sec ond comparison project, including gas dynamics ( llliev et al] 
120091) . 

Adding helium implies introducing another source of 
frequency-dependent opacity, thus making a multi-frequency ap- 
proach inevitable. In addition the on-the-spot approximation be- 
comes more complicated as one has to take into account how re- 
combination photons from helium affect the hydrogen ionization. 
As one moves to higher photon energies, one should also take into 
account the secondary ionizations caused by the superthermal elec- 
trons produced when high energy photons ionize the atoms and 
ions. Including EUV and soft X-ray (SX) photons therefore is a 
non-trivial extension of the photo-ionization calculations which we 
describe and test in this paper. 

In terms of the physical processes included the method we 
present here is similar to a number of others published in re- 
cent years, but diffe rs in the algorith ms used. The cod es CRASH 
dMaselli et al.l2009h and LICORICE teaek et alj201f3) use Monte 
Carlo techni ques. The codes SPHR AY dAltav et alT 120081) and 
TRAPHIC jPawlik& Schavd l201ll) implement ray tracing on 
particle data whereas RADAMESH JCantalupo & PorcianfeOllh 
uses an adaptive mesh approach. 

The lay-out of the paper is as follows. In Section [2] we give 



an overview of the basic algorithmic ideas behind C 2 -RAY to then 
proceed in Section[3]with an overview of how we extended its capa- 
bilities to handle harder photons. Section|4]contains the description 
of a series of one and three-dimensional tests for the new method, 
also evaluating the effects of the various new elements such as sec- 
ondary ionizations and the coupled on-the-spot approximation. A 
series of appendices describe several important elements in more 
detail. 



2 REMINDER OF BASIC STEPS OF THE ORIGINAL 
C 2 -RAY ALGORITHM 

The C 2 -RAY method was developed to be a time-dependent photo- 
ionization algorithm that could be efficiently combined with a 
hydrodynamics calculation, and not impose impractically short 
timesteps and small cell sizes on the latter. This is achieved by 
assuming that the ionization evolution of individual cells follows 
an exponential decay to the equilibrium solution and that a time- 
averaged value of the optical depth can be used to describe the ef- 
fect of a cell on the radiative transfer during the entire timestep. 
This approach is able to correctly track the progress of ionization 
fronts over many cells during one timestep. In addition, optically 
thick cells are dealt with by defining the photo-ionization rate such 
that it is consistent with the number of photons absorbed inside a 
cell. C 2 -RAY was described and tested in detail in M+06. Here we 
summarize some of the ideas in order to define our notation and 
provide an introduction to the extensions described in Sect. [3] 

The evolution of the ionized hydrogen fraction derives from 
the set of chemical evolution equations: 
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where xm is the neutral hydrogen fraction, xuu is the ionized hy- 
drogen fraction, n e is the electron density, Thi is the hydrogen 
photo-ionization rate (see below), Chi is the collisional ionization 
rate and a^n is the recombination rate. Since xm + xmi = 1, we 
obtain 

-t-^hii = — (Thi + C'mn e + a^ u n e ) xmi + (Tm + C'mn e ) ■ 
at v .. ' v .. - 



A H 



(2) 

Here we introduce the notation, A H and <? H , which allows us to 
write the equation in the general vector form, useful later on, 

d 



— x = A x + q . 
dt a 



(3) 



As is well known, the general solution x(t) to a set of equations of 
this type is the sum of the solution to the homogeneous case, Xh{t), 
where g — o, and a particular solution x p : 



x(t) — Xh(t) + x p with Xh = Cj Xj e* 



(4) 



Here, n is the rank of A, i.e. the number of coupled equations (1 
in the case of hydrogen), Ai are the eigenvalues of A, X; are the 
corresponding eigenvectors and c; are coefficients which can be 
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calculated from the boundary condition x(t = 0) = Xo: 

71 



£C () 



(5) 



In subsequent timesteps, xo is the state at the end of the previous 
timestep. In the case of a constant g, the particular solution can be 
the equilibrium solution given by 

Ax p + g — o . 

For the simplest, hydrogen only case, the coefficients are thus 
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Here, we added the superscript H for hydrogen and we skipped the 
subscript 1 since for hydrogen, n = 1. 

The photon-conserving photo-ionization rate V in each cell 
used in A is calculated using 
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where (At„} is the time averaged optical depth over the cell and 
V^heii is the volume of the shell the cell belongs to. This quantity 
can be calculated from the time evolution of the neutral fraction 
(Eq.[4]l. By solving these two equations l[4]and[7]l in alternating or- 
der one iterates to convergence as illustrated in Fig.[T] This iteration 
also involves the electron density n e which is calculated from the 
time averaged ionized fraction. The time averaged optical depth to 
the cell (t v ) is calculated by short characteristic ray-tracing over 
the solutions found for cells lying nearer to the source. This makes 
the algorithm causal. 

In the case of multiple sources, the iteration as shown in Fig.Q] 
is split up in two parts: The first part, including step three (find- 
ing the ionization and heating rates in each cell) is done for each 
source, looping through the entire computational grid using short 
characteristic ray tracing. For each cell, the rates from all sources 
are added. These total rates are used in the remaining two steps in 
the iteration. See also M+06 for a flow chart and a description of 
the implementation how to loop through the source list. 

Note however that the flow chart in M+06 for the single source 
loop (figure 4) incorrectly includes the last two steps of the single 
cell loop. The electron density and the time averaged ionization 
fractions are in fact not updated in the source loop but are updated 
first after the photo-ionization rates from all sources are summed to 
a global photo-ionization rate. 



3 EXTENDING C -RAY 

The original C 2 -RAY methodology works well for soft, stellar 
spectra. In this case, the IGM can be considered to be hydrogen 
only since there are not many photons capable of ionizing He II 
and helium can be assumed to be singly ionized everywhere where 
hydrogen is ionized. In Sect. 14.3.11 we show that the morphology 
of the ionization fraction field in a cosmological reionization sim- 
ulation with stellar sources (only), is indeed hardly affected by the 
inclusion of helium. However, when the spectrum has a significant 
amount of SX photons, helium contributes significantly to the op- 
tical depth and a multi-frequency approach is required: at frequen- 



Find column density to current cell 



Find optical depth over current cell 



Find photo-ionization and heating rates 
in current cell (Eq. 7) 



Find electron density in current cell 



Find new ionization fraction 
(+time averaged ) in current cell (Eq.4) 




Figure 1. Iteration scheme of C 2 -RAY for a single cell, conceptionally as 
in M+06, figure 2 



cies higher than the ionization threshold of Hell, the ionization 
cross-sections of Hel and Hell are roughly an order of magnitude 
larger than the HI ionization cross-section. Therefore, neglecting 
helium in the case of sources with hard spectra will underestimate 
the optical depth substantially. We therefore have to add helium 
chemistry to C 2 -RAY. 

In order to include helium, both the chemical evolution equa- 
tion, Eq. ©, and the calculation of the ionization rate, Eq. ((7), have 
to be changed. Additionally, as shown below, the iteration scheme 
from Fig.Q]has to be modified. We describe each of these changes 
here. 



3.1 Chemical evolution equation 

The procedure for adding helium photo-ionization to our calcula- 
tions is by itself relatively straightforward as the basic algorithmic 
idea described in Sect. [2] provides the frame work for this. How- 
ever, a complicating factor is the presence of ionizing recombina- 
tion photons since they couple the rate equations of hydrogen and 
helium. Here we present two approaches for dealing with these, 
where the first one is less accurate, but simpler and more similar to 
what other authors have used. We compare the results of these two 
approaches in Sect. [4] 



3.1.1 Simple recombination: no coupling of species 

When ions recombine, photons are emitted. In case of recombina- 
tion of hydrogen ions, only recombinations to the ground state re- 
sult in photons energetic enough to ionize hydrogen. If one assumes 
these photons to ionize immediately another hydrogen atom close 
by, this i s called the on the spot (OT S) approximation for hydro- 
gen (e.g. lOsterbrock & Ferl and 2006). In this approximation, the 
recombination coefficient to all states of hydrogen, ct A is replaced 
by the recombination coefficient to all states but the ground state, 
a B . For a mix of hydrogen and helium, the OTS approximation is 
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more complicated as helium recombination photons can ionize both 
hydrogen and helium. However, as a first step, we assume that pho- 
tons from recombinations to the ground state can only ionize the 
same species from which they originate and that in recombinations 
to other states than the ground state no ionizing photons are emit- 
ted. That means, we use the a B recombination-coefficients for all 
species. In the following we refer to this as the "uncoupled on-the- 
spot approximation" (U-OTS). In this approximation hydrogen and 
helium can be treated separately. For helium, the set of chemical 
evolution equations (in analogy to Eq.Q]for hydrogen) is: 



di 




n e O!HcII 
-n e Q!HcII — £/HcII W e aHeIII 
t^Hell — ^eQHcIII 



where xuai is the neutral helium fraction, XHcii is the singly ion- 
ized helium fraction and £hcIii is the doubly ionized helium frac- 
tion. We grouped the ionizing 'up rates' into one term 
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The subscripts on C, a and F indicate on which species they act, so 

a HeIII 

He II ^ He III. Using the fact that xna + xmn + ^Hciii = 1, 

r Hen 

the equivalent equation to Eq. l[2j can be written as 
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where the vector g and the matrix A have the following forms 
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(11) 



The solution of this set of linear differentia l equations for the U- 
OTS case can be found in lAltav et alj |2008) or with the here intro- 
duced notation in Appendix [A] 



3.1.2 On the spot approximation: coupling of species 

In reality, recombination photons from helium can ionize either hy- 
drogen or helium, introducing the need to couple the rate equations 
for the two elements. This is the proper OTS approximation. Ta- 
ble Q] gives an overview of the recombination processes affecting 
hydrogen and helium fractions in the OTS approximation. 

To implement the OTS approximation we follow 
lOsterbrock & Ferlandl d2006h for dealing w i th the recombina- 
tions of Hell to Hel and F lower & Perinottol Jl980h for those of 
Helll to Hell (except for recombinations to the ground-state). 
Mostly we also use their notation. 

For hydrogen, we take the a s -recombination coefficient a ml . 
For Hell, the photons from recombinations to the ground state are 
distributed between helium and hydrogen depending on the frac- 
tion of optical depth at the helium ionization threshold frequency, 
fjh • a fraction y goes into hydrogen ionization, a fraction 1 — y 
goes into helium ionization. The photons from states other than the 
ground state contribute with a fraction p to hydrogen ionization. 
Similarly, for Helll the recombinations to the ground state ionize 



Hell, Hel and HI, depending on the relative optical depth over the 
cell in question at i/jJJ . Of those recombination photons, a fraction 
y% goes into Hell ionization, a fraction y\ g° es into Hel ionization 
and a fraction 1 — y^ — y\ goes into HI ionization. Here, the frac- 
tions y, y%, and y\ are dependent on the relative optical depth of the 
species at the threshold frequencies of Hel and Hell, as described 
below. 

Recombinations of Helll to other states than the ground state 
contribute to ionization of HI and Hel. Those recombinations lead 
either to two photon emissi on (in a fraction v of the cases, where 
v is temperature dependent. Hummer & Seatonl 19641) . of which on 
average a fraction I is energetic enough to ionize HI and a frac- 
tion m is energetic enough to ionize Hel; therefore, a fraction 
vw = v (I — m + my) goes into HI ionization and a fraction 
v (m (1 — y)) goes into Hel ionization. The remaining fraction, 
1 — v, leads to emission of a He Lyman a photon. Those photons 
are absorbed by any species to a fraction /. By letting escape some 
of the helium Ly a photons (/ 7^ 1) we lose them since we do 
not include those at larger distances from the source. Of the ab- 
sorbed He Lyman a photons, a fraction z goes into HI ionization 
and the remaining fraction 1 — z into Hel. Additionally, the Balmer 
continuum emission photons (QhcIii) can ionize hydrogen. TableQ] 
summarizes the photon emitting recombination processes included 
in the on the spot treatment. The numerical parameters used are 
listed in Table[2] 

We can now introduce these terms in the general equation, 
Eq. {3}: 
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The density fractions nHc/«H in A12 and A13 are due to the 
fact that the equations are written in terms of fractions, not in terms 
of densities. The complete solution for this set of equations is pre- 
sented in Appendix |B1 

3.2 Extending the iteration mechanism 

In Sect. 13.1.21 we treat the set of chemical evolution equations 
dx/dt — f(x) as a set of linear differential equations, Eq. (|3}, 
and use iterations to obtain the correct electron density and ioniza- 
tion rates. However, further dependences on the different ionization 
fractions are hidden in the parameters y, z, w, y% and y\. We found 
that this hidden extra non-linearity can complicate the convergence 
of the iteration. We therefore extended our iteration scheme to make 
it more robust. The extension consists of updating the parameters 
y, z, w, j/2 an d J/2 a f ter the chemical evolution equation has been 
solved and then solving the chemical evolution equation for a sec- 
ond time with this new set of parameters. We then take the mean 
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Table 1. Summary of recombination processes included in the OTS treatment 
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Table 2. Overview of the numerical parameters used in the OTS approximation 



p = 0.96 or 0.66 depending on n e JOsterbrock & Ferlanol2 006). in our cases of interest (cosmological simulations) always 0.96 

y = tjj/(tjj + THoi) at i/^ cI (ionization threshold of He I) 

l/f = T Hell/( r H + T Hel + THell) at ^th (ionization threshold of He II) 

v\ = T Hcl/( r H + T Hcl + T Hcll) at f^ 11 (ionization threshold of He II) 

z = th/(th +t Hc i) at hi/ = 40.8eV (Hel Ly a) 

f = 1 to 0. 1 ("escape" fraction of Ly a photons) depending on the neutral fraction 

v = temperature dependent coefficient! Hummer & Seaton 1964) 

ui = (I — m) + my 

I = 1.425, fraction of photons from 2-photon decay, energetic enough to ionize hydrogen jFlower & Perinottdll980l) 

m = 0.737, fraction of photons from 2-photon decay, energetic enough to ionize neutral helium (Flower & P erinot to| ll980l) 



of these first and second solutions for the time-averaged ionization 
fractions and use this mean to calculate the photo-ionization rates. 
These extra steps mean that the iteration scheme now consists of 
seven instead of four steps, see Fig. [2] 



3.3 Calculating the photon rates 

For a pure hydrogen medium Eq. (|7) gives the photo-ionization 
rate, i.e. F = Fhi and r = Tin. In the case of a medium consisting 
of hydrogen and helium, this simple treatment can only be used for 
photons with frequencies below the ionization threshold for neu- 
tral helium, v^ 1 . Above the ionization threshold frequency of ion- 
ized helium, all three species, HI, Hel and Helll contribute 
to the optical depth. In the frequency bin between those threshold 
frequencies, HI atoms and Hel atoms contribute to the total op- 
tical depth. Therefore, the minimum number of frequency bins to 
consider separately when calculating the photo ionization rates for 
each species is three, henceforth referred to as (frequency) bins. As 
we explain below, we use a single frequency dependence for all 
species in each bin. Since the the different frequency dependences 
of the ionization cross sections a of the species are very different 



in the bins, it is useful to further subdivide bin 2 and bin 3, hence- 
forth referred to as (frequency) sub-bins. To refer to a particular 
sub-bin, we use the following notation: n.m refers to sub-bin m of 
bin n (n = 2, 3). In Ap pendix ICl we show o ur power-law fits to the 
cross-section data from Ver ner et al. (1996). 

The general treatment in every such sub-bin follows M+06 
and is schematically illustrated in Fig. [3] To avoid expensive inte- 
grations for each photo-ionization calculation, we tabulate the total 
ionization rate r to t as a function of total optical depth r to t at the 
minimum frequency of the sub-bin in question, assuming for all rel- 
evant species the same frequency dependence for the cross-se ction, 
following the concept outlined in lTenorio-Tagle et al.l a 1985b . The 
assumption of having the same frequency dependence allows the 
summation of the optical depth of each species and the tabulation 
of the photo-ionization rate as a function of this single total optical 
depth. Since the frequency dependences of the cross-sections of the 
different species are in fact not identical (see for example Fig. lC 1 1 >. 
this is an approximation. The more sub-bins we use the more ac- 
curately we follow the actual shape of the (frequency dependence) 
curves of the different cross-sections. The frequency dependence 
imposed on all species is that of cthcI for all sub-bins of bin 2 and 
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Find total column density to current cell 



Find optical depths over current cell 



Find photo-ionization and heating rates 
in current cell (Eq. 7/13) 



Find electron density in current cell 



Find new ionization fractions 
(+time averaged ) in current cell (Eq.4) 



z 



Find optical depths over current cell to 
update y,z,w, y2a and y2b parameters in 
current cell 



Find new ionization fraction 
(+time averaged ) in current cell (Eq.4) 



Average from ionization fractions 
calculated in Step A and B 




Figure 2. Extended iteration scheme for C 2 -RAY including helium with 
coupling of the species. 



that of (THeii for all sub-bins of bin 3. The table entries are the in- 
tegrals of ionization rate over the frequency range of the sub-bin in 
question, thus the total ionization rate due to photons with frequen- 
cies in that frequency sub-bin. 

In order to use this pre-calculated table and determine the total 
photo-ionization rate in a given cell, we need to determine the ion- 
ization cross-sections for every species at the minimu m frequency 
of eve ry sub-bin. We use the fitting formula from IVerner et alj 
dl996h for the cross-sections and compute the total optical depth 
Ttot as the sum of the optical depths of all species at the minimum 
frequency of each sub-bin. This optical depth is used to determine 
the total photo ionization rate in the table. 

Next, this total photo-ionization rate has to be split up into 
ionization rates for HI (Thi), Hel (Thci) and in case of bin 3, Hell 
(THeii). This is done according to the relative fractions of optical 
depth in each frequency interval. 



, Ttot = > ' 



Ttot 



HI, Hel and Hell (14) 



In Appendix |D] we evaluate other approaches for this distribution 
that have been proposed in the literature. 

We calculate these fractions of optical depth at the minimum 
frequency of each sub-bin. This implies that we also here assume a 
single power law fit for the frequency dependence of all species in 



each sub-bin. Doing so overestimates the contribution of hydrogen 
to the optical depth in frequency bin 2 since the HI ionization cross 
section drops faster with frequency than the ionization cross sec- 
tion of Hel. This leads therefore to an overestimate of Thi/ThcI- 
In frequency bin 3, the slopes of the cross-section curves are in gen- 
eral more similar (see Appendix [Cj but still steeper for HI which 
results in an overestimate of Thi/ThcH- The more sub-bins used, 
the smaller the overestimates will be. We test the convergence of 
this in Appendix [E] We find that increasing the number of sub-bins 
in bin 2 improves the slope of the ionization front while increasing 
the number of sub-bins in bin 3 improves the ionization fractions 
inside and outside the front. 



3.4 Heating and secondary ionizations 

One of the motivations for extending C 2 -RAY, is to investigate the 
effects that helium and hard spectra have on the temperature evo- 
lution of the IGM during the EoR. In this section we first describe 
the implementation of the temperature calculation. This is followed 
by remarks on the inclusion of secondary ionizations and consider- 
ations about the additional timestep restrictions caused by the tem- 
perature calculation (in addition to the ionization calculation). 

We use the ideal gas law PV = NksT, where ks is Boltz- 
mann constant, to calculate the pressure P from the temperature 
T and the gas particle number density N/V — nu + TiHe + n e 
in each cell of volume V. From the pressure we then calculate 
the internal energy per volume V (in each cell), using the di- 
mensionless heat capacity (per particle) at constant volume, cv' 
u int _ jytntyy _ Cy p p or a monatolruc gaSi Cy = 3/2. So the 

equation to convert temperature into internal energy density reads: 

int. 3 . 



-fesT(nH + nHc + n e 



(15) 



This energy density is affected by the heating (H) and cooling (C) 
of the gas per (cell-) volume and per unit time. In order to follow 
the temperature evolution of the gas, we therefore have to solve the 
equation 



dt 



= H-C. 



(16) 



As the only contribution to the heating rate W, we consider 
photo-ionization heating. In every photo-ionization of species i, the 
excess energy, (v — vth(i)) is transferred to the electron released. If 
one assumes that the cells are optically thick, (almost) all photons 
are absorbed and the average energy per photo ionization is simply 
the average excess energy over the whole spectrum. This was the 
approach used by the original C 2 -RAY in the tests with temperature 
evolution in 1+060 

To properly take into account the effects of different optical 
depth on the heating we can calculate "H in analogy to the ioniza- 
tion rate V (Eq.[7). It should be remembered that, although we write 
this as one equation, for the photo-ionization rate it is in fact the 
difference between the ingoing photo-ionization rate (calculated on 
the basis of (t)) and the outgoing photo-ionization rate (calculated 
on the basis of the (t) + (At)). For the heating, these quantities 
are rather abstract since they symbolize the excess energy still in 
the form of photons when entering the cell and the amount of ex- 
cess energy still in the form of photons leaving the cell. Neverthe- 



1 The H-only C 2 RAY used in iMellema et ai] l2006ah and lArthur et al] 
feOl ll) actually used a one species, one frequency bin version of the heating 
method described below. 
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STEP 1: calculate t as the sum of t , t and t at v . for every sub-bin 
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STEP 2: look up T(t) 
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STEP 3: in each sub-bin 
'„,., n calculate f based 

> f 3 ° on the realtive t(v . ) 



f 1 x r 



f 3 1 x r 



f 311 x r 



. r . 3 ° 



tot STEP 4: for each species: 
Add the contributions to T 
from all the sub-bins 



HI 



Figure 3. Sketch showing how to calculate the ionization rate contribution for bin 3. In step 1, the total optical depth in each sub-bin is calculated. In step 

2, the pre-calculated photo ionization rate table is used to calculate the total photo ionization rate in each sub-bin. On the basis of the optical depth of each 
species at the minimum frequency u mln in each sub-bin, the fractions fi going into hydrogen-, helium- and ionized helium- ionization are calculated in step 

3. Finally, the products of the fractions and total photo ionization rates are summed over all sub-bins to result in the ionization rate for each species, Th, ^Hel 
and r Hc ll- 



less, what is tabulated as function of optical depth is this excess 
energy as a function of optical depth. The corresponding equation 
to Eq. {7} is then 



H(i) = 



h{v - f th (i))L„e- <T - > 1 - e -< AT "> 



y(sub — bin) 



hv 



du 



(17) 



for species i. Since the excess energy is different for each species, 
instead of having one table of heating rates in each sub-bin, we 
need three tables, where the excess energy is with respect to the 
threshold frequency of ^ t h(HI), i/ t h(HeI) and ^ t h(HeII). The to- 
tal heating in each frequency sub-bin is naturally analogous to the 
photo-ionization given by 



n = y —u(i) . 



(18) 



Some of the electrons that received excess energy (y ~ ZAh(*)) 
after being released from the bound state of atom i, will collide 
with bound electrons of other atoms (or ions) rather than other free 
electrons. This transfers energy to the bound electron. If the trans- 



ferred energy is greater than the binding energy of this electron, 
this electron is released. This process is called secondary ioniza- 
tion. The probability of such an event depends on the energy of 
the primary electron and on the ionizatio n state of the gas. The 
full pr o cess requires carefu l modelling, but Shull & va n Steenberd 
( 1985). lRicotti et al.l ( T20o3) and lValdes & Ferrar ij fc008h . to name 
a few, published separable functional relations dependent on pri- 
mary electron energy and hydrogen ionization fraction to describe 
what fraction of the energy deposited in primary electrons goes 
into secondary ionizations of either HI or Hel, /™ HeI , and what 
into heat, /h ca t- In these relations it is assumed that xnu = ^Heii 
and secondary ionizations of Hell are n eglected. We i mplem ent the 
separable functional relationship as in iRicotti et alj J2002T) which 
converges for high electron energies to the functional form of 
IShull & van Steenberd ( fl985T) . 

iFurlanetto & Stoevetl bold) pointed out that there is in gen- 
eral no simple separable functional form for deciding the fractions 
that go into secondary ionizations and heating. Rather the depen- 
dence of the relative fractions varies with ionization fraction in a 
complex way, see their figure 5. Ho wever, the complex relation 
given by IFurlanetto & Stoevetl feOlfj) still assumes shii = ZHeii 
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and no secondary ionizations of Hell. Given these limitations we 
decided that at this point there is no clear benefit in implementing 
these more computationally expensive relations. 

For the cooling rate, we include free-free and recombination 
cooling for H II, He II and He III and collisional excitation cooling 
for H I, He I and He II. A full overview of the radiative cooling 
rates used is given in Appendix [G] In case of cosmological simu- 
lations, cosmological cooling due to the expansion of the universe 
and Compton cooling against the cosmic microwave background 
photons are included as well. 

In order to numerically solve Eq. d!6l > we use forward Eu- 
ler integration. However, the cooling rate depends sensitively on 
the gas temperature which is changing because f the heating. To 
accurately follow this behaviour we use the forward Euler method 
with sub-timesteps determined by a limit on the temperature change 
(typically 10%). We keep the heating term H (which represents the 
average heating rate over the timestep) constant but use the temper- 
ature from the previous sub-timestep to calculate the cooling rate 
C. 

Above, we described how to calculate the heating rate in anal- 
ogy to the photo-ionization rate. Specifically, we use the time- 
averaged optical depth (to the cell and in the cell) both in Eq. (7} 
and Eq. dl7| (. There is however an important difference between the 
dependence of the photo-ionization rates and the heating rates on 
optical depth. Although Eq. 10 gives the correct rate of absorbed 
photons, for the heating the frequency of the photons matters. Let 
us first consider the source cell, i.e., the optical depth to the cell,r, 
is constant with time and equal 0. While for the photo ionization, 

r«r) At ) = ^ — ^ " , (19) 

for the heating in general 

f^HMt'^dt' 
H«r) At ) < 4> y " , (20) 

because of the (y — fth(i)) term in the integral for T-i. Qualitatively 
one could say that in the beginning of the timestep the heating per 
photo-ionization is the optically thick value (average photon energy 
over the whole spectrum), for later times it shifts to lower values, 
the limit of which is given by the optically thin value (average pho- 
ton energy over the spectrum weighted with the photo-ionization 
cross-section) 

For this first cell where the only dependence is on the optical 
depth inside the cell, it would be possible to tabulate a correction 
factor based on the change of the optical depth during a timestep. 
However, in the more general case where the optical depth to the 
cell depends on the time varying optical depth of the other cells 
along the ray, a local fix (such as using an energy- average based 
on the change of optical depth in the cell or sub-timestepping on 
a cell-by-cell basis) to the heating rate is no longer possible. We 
explored different approaches for calculating accurate heating rates 
using large timesteps, but we were unable to find a general solution. 

In Appendix|F]we investigate how the heating depends on the 
choice of the timestep. From these tests we conclude the following: 

• If we are interested in time scales larger than the recombina- 
tion time, we can still use large timesteps, as the initial heating is 
no longer dominating. 

• If the cooling time t coo \ < At we can also use large timesteps, 
as the temperature is set by the equilibrium heating and cooling 
rates. This is the case for typical interstellar medium conditions. 

• If we are interested in time scales below the recombination 



Table 3. Parameters for testing the code. 



tih /cm 


constant hydrogen number density 


n H Jcm~ 3 


constant helium number density 


N-y/S 


rate of photons in the energy interval 




[13.6,5441.6] eV 


T cff /K 


effective temperature of the black body source in 




case of BB source 


P 


power law index in case of PL source 


T ini ,x 


initial temperature and initial ionization states, H II, 




He H and He III 


Ar 


the cell size 


At 


the timestep 




number of sub-bins in frequency bins 2 and 3 


OTS/U-OTS/ a A 


assuming the (U-)OTS-approximation, or using 



a A recombination rates 



time, an accurate value for the temperature requires timesteps of the 
order of the ionization time. This constraint becomes more strict if 
the cells are very optically thick. 

Given these conclusions it is difficult to provide a simple 
recipe for choosing the timestep. We therefore recommend testing 
for numerical convergence if accurate temperatures are required. 



4 TESTING THE CODE 

To test the validity of the approximations we made for the sake of 
code efficiency, we performed a series of tests. First of all we vali- 
dated the new version against the old version of C 2 -RAY by setting 
the helium abundance to a very low value. This test showed that 
the new version has the same photon-conserving properties as the 
method presented in M+06. For testing the time dependent solution 
with helium, including the OTS approximation, secondary ioniza- 
tion and temperature evolution, we would need a fully validated 
time dependent photo-ionization code, which as far as we know is 
not publicly available . W e were therefore forced, just as for exam- 
ple lBaek et al] d2010t> and lPawlik & Schavel ( 120111) to compare our 
results against the one-dimensional photo-ionization equilibrium 
code CLOUDY (version 08.00, last described in Ferla ncfet alj 
1998). Below we present the results of two sets of such one- 
dimensional tests, one set in which the source has a black body 
(BB) spectrum with an effective temperature of T e s, the other for 
a source with power-law (PL) like spectrum, where the energy- 
distribution can be described as L(v) tx . For the BB source 
we test various aspects of the calculation of the ionization fractions 
while keeping the temperature constant. For the PL case we also 
consider the temperature evolution. We chose this approach since 
the latter case having a wide energy range of hard photons, consti- 
tutes a more difficult test for the photo-heating. Table [3] gives an 
overview of the parameters used for all the one-dimensional tests. 
To demonstrate the multi-dimensional and multi-source capabili- 
ties of the extended C 2 -RAY we also present the results of two test 
problems using three-dimensional cosmological density fields. 

4.1 Test-suite 1: Black body source 

The first set of tests considers the expansion of an ionized region 
produced by a single source into a constant density medium with 
constant temperature. Under the assumption of spherical symmetry 
this can be calculated in one dimension for the radial distance r to 
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the source. The parameters that are not varied in this set of tests are: 
n H /cm" 3 = 1.0 x 10~ 3 , n He /cm" 3 = 8.70 x 10~ 5 , 7V 7 /s = 5.0 x 
10 48 , T off /K = 10 s . T ini = 10000K, x = (10" 4 °, 10~ 40 , 10" 40 ) 
(i.e. completely neutral), Ar/pc= 150, Af/yr=10 7 and (712, n^) — 
(10,14). 

We choose these parameters so that the helium fraction fn e = 
TiHe/(^H + ?iHc) = 0.08 and the remaining physical parameters 
are the same as in Test 2 of 1+06, except that we do not follow 
the temperature evolution here. The timestep At is ~ 10% of the 
recombination time in a fully ionized medium, tihQ; A ' b , and ~ 10 4 
times the ionization time for the first cell. 



4.1.1 Test 1 A 

In the first test of this suite we assume that all ionizing photons from 
recombinations can escape, which means that we are using the a A 
rates. This test is meant to show how well our basic approach can 
match the equilibrium solution from CLOUDY. In Fig.|4]we show 
the results together with that equilibrium solution. It can be seen 
that after t = 10 10 yr the general agreement with the CLOUDY 
equilibrium solution is excellent for distances below 12 kpc, but not 
for larger distances. However, since the recombination timescale is 
given by (n e a)~ , it follows that for ionization fractions below 
~ 1%, the recombination timescale is larger than ~ 10 10 yr. In 
other words, these outer regions have not yet reached equilibrium. 
Another noticeable feature in Fig.|4]is that the Helll fraction shows 
a transient bump around r = 4 kpc. This bump disappears as the 
equilibrium solution is approached, but both, our results and the 
equilibrium solution still show a slope change in the Helll curve 
around that same distance. From all this we conclude that the re- 
sults of this test show that our basic multi-frequency radiative trans- 
fer is consistent with CLOUDY's. 



4.1.2 Test IB 

In the second test, TEST1 B, we use the on the spot approxima- 
tion as described in Section [3. 1.21 and otherwise set the same pa- 
rameters as in the previous test. We show the ionization profiles 
in Fig. [5] together with two results from CLOUDY, one using the 
OTS approximation (green dashed line) and the other using full 
radiative transfer of recombination photons (blue dashed line). As 
is apparent in Fig. [5] the agreement between the C 2 -RAY and the 
OTS solution of CLOUDY is good for those distances where the 
equilibrium solution has been reached (r < 12 kpc at t — 10 10 yr, 
just as in TEST1 A). The general pattern of evolution of the ion- 
ization fractions is quite similar to the one in TEST 1 A, although 
the OTS approximation does of course change the detailed values 
of the fractions. 

The CLOUDY solution using full radiative transfer of recom- 
bination photons can be used to evaluate the validity of the OTS 
approximation for this particular test problem. Inspection of Fig. [5] 
shows that the error introduced by using on the spot approxima- 
tion is barely noticeable, but larger for hydrogen and largest for 
the HI fraction inside the HII region. Closer inspection shows that 
the neutral hydrogen profile inside the HII region, up to a neutral 
fraction of about 1%, follows the curve from the a A recombination 
coefficient from TEST1 A. This is not surprising since most of the 
recombination photons from this highly ionized region will not be 
absorbed on the spot, but escape to larger distances. 

To gauge the importance of applying the coupled OTS rather 
than the popular U-OTS (i.e. using a B rates, see Sect. 13. Lit , we 



also plot the ionization profiles after 10 yr for the latter ap- 
proach (thin black lines). It can be seen that the differences be- 
tween the OTS and U-OTS results are larger than between the OTS 
approximation and the full radiative transfer results, when compar- 
ing the (almost) equilibrium solutions. From this we conclude that 
although more complicated to implement, the coupled OTS approx- 
imation should be the preferred approach. 

Since C 2 -RAY cannot (now) deal with the diffuse photons in 
another way than using the OTS approximation, we are not able 
to investigate the effects of using this approximation d uring the 
growth of the H II region. Cantalupo & Po rcianil d201ll) compare 
results from what we call the U-OTS approximation with full ra- 
diative transfer of diffuse photons using their code RADAMESH. 
They find rather large effects at times when the equilibrium solu- 
tion has not yet been reached. This warrants further investigation, 
probably best pursued in manner similar to 1+06, involving results 
from multiple codes. 

4.2 Test-suite 2: A power law source 

In a second set of tests, we use UV-photon emitting sources with 
a power-law spectrum, L(y) oc implying that the number of 
photons goes as f(y) oc v~^ +1 >. For these we considered the two 
cases B — 1 and B = 2. We present the results for simulations with 
temperature evolution using Ti n i = 10 2 K. We only use the full 
OTS approximation, apply the secondary ionizations and choose 
the number of sub-bins in bin 2 and 3 to be (712, n.3) = (26, 20). 
The remaining parameters are as in TEST 1. 

We found a timestep of At = 10 yr to be sufficient to obtain 
convergence in the temperature evolution. Choosing a timestep of 
the order of the ionization time of the first cell (At — 10 3 yr) 
gave temperature results which for cells close to the source were 
only different by at most 7%. For the parameters of this test, the 
optical depth at the ionization threshold for hydrogen for one cell 
is r ~ 3.2, in between the two cases tested in Appendix [F] and so 
this result for the timestep is consistent with those. 

We show the profiles for HI/HII and Hel/Hell/Helll and T 
for three times and the two different power law spectra in Fig. [6] 
We also plot the equilibrium solutions from CLOUDY with OTS 
for those two spectra. We find that at t = 10 9 yr the C 2 -RAY re- 
sults are close to the CLOUDY ones for those distances where 
the equilibrium solution has been reached, although not as close 
as for the results of TEST 1. These somewhat larger differences 
with the CLOUDY results cannot be explained by effects from the 
temperature calculation since we found comparable discrepancies 
when imposing a constant temperature. However, overall the match 
is still reasonable. 

The time evolution shows that the two values of B produce 
similar results, with initially quite steep profiles for the ionization 
fractions for all times and a relatively steep front in the temperature 
evolution even for the quite hard B = 1 spectrum at the relatively 
late time of t = 10 7 yr. Comparing the B = 1 and B — 2 results 
it can be noticed that the former gives a higher degree of hydrogen 
and helium ionization outside the front than the latter. Further it can 
be seen that the residual neutral hydrogen fraction inside the front 
is higher for the results with /3 = 1. Just as in TEST1, the bumps 
in the Helll fraction are transient phenomena, although even the 
equilibrium solution shows a slope change around a distance of 
~ 6 kpc. 

iTittlev & Mei ksin (2007) reported the HII front to be trailing 
behind the HcIII front for spectra with B < 1.8. However, even in 
our B = 1 results we find the HII fraction to be always larger than 
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Figure 4. Results from TEST1 A. Fractions of HI (upper left panel), HII (upper middle panel), Hel (lower left panel), Hell (lower middle panel) and Helll 
(lower right panel) at times t/yr = [1 X 10 7 , 1 X 10 8 , 1 X 10 9 , 1 X 10 10 ] as indicated by line thickness in the legend. We also show the equilibrium solution of 
CLOUDY (green dashed). Note that due to the low electron density in the outer regions, the equilibrium solution has not yet been reached beyond a distance 
of about 12 kpc at t = 10 10 yr. 



the Helll fraction. We were able to obtain (double) front crossings 
when using the a A recombination rates (no OTS) and disabling 
the secondary ionizations. However, the locations where the Helll 
fraction was found to be larger than the HII fraction were in a lim- 
ited area in the ionization fraction-time space, roughly in the in- 
terval [0.1, 10 8 yr] to [0.01, 10 9 yr]. We therefore conclude that 
trailing HII fronts are at best a marginal effect. 

In order to evaluate the effects of secondary ionizations we 
also ran TEST 2 without them. We found that this led to changes 
larger than the differences we found between our results and the 
equilibrium solution of CLOUDY. We therefore conclude that it is 
worth to include secondary ionizations in the calculation. 



4.3 Cosmological tests 

In this section, we test the 3D-version of the extended C 2 -RAY on 
two cosmological density fields. The first test is a larger scale cos- 
mological test without temperature evolution to evaluate the effect 
of helium at a fixed temperature in a simulation with sources with 
a soft spectrum. The second test is a rather small cosmological vol- 
ume but includes temperature evolution: we redo the cosmological 
test-problem with multiple sources from 1+06 (their test4) to test 
the effect of helium on the heating and on the hydrogen ionization 
fraction field. 



4.3. 1 TEST 3: Effect of helium on the morphology of the 

hydrogen ionization fraction field during EoR without 
temperature evolution for stellar type sources 

Many cosmological reionization simulations only include hydro- 
gen and implicitly assume that helium is singly ionized every- 
where where hydrogen is ionized. The used hydrogen number den- 
sity is therefore equal to the total number density. In this sec- 
tion, we test if the morphology of H II regions in a reioniza- 
tion simulation with stellar sources changes if helium is included. 
For this compa r ison, we u se simulati o n 53 Mpc_g8,7_130S from 
iFriedrich et al J d201lh and llliev et al] d201ll) . The electron scat- 
tering optical depth produced by this simulation, r cs = 0.083, is 
consistent with the 1-cr range allowed by the seven year WMAP 
results, t cs = 0.088 ± 0.015 dKomatsu et~al]l201 ll) . In this simu- 
lation, the number of ionizing photons produced by a dark matter 
halo of mass M is defined through 

N -y = 9j , nn - , (2i) 

lus l m m p 

where N y is the number of ionizing photons emitted per Myr, 
fib = 0.044, fi m = 0.27 and m v is the proton mas s (This equa- 
tion included incorrectly a /j, in lFriedrich et al Jboill . equation 1). 
Massive halos are assigned an efficiency g 7 = 8.7 while low mass 
sources have an efficiency of g-, = 130 and are suppressed in re- 
gions where the ionization fraction (of hydrogen) is higher than 
10%. To evaluate the effect of helium on the morphology of H II 
regions we use the dimensionless power spectrum of the H II frac- 
tion A*a In Fig. [7] we show the power spectra for this simulation 
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Figure 5. Results from TEST1 B. Fractions of HI (upper left panel), HII (upper middle panel), Hel (lower left panel), Hell (lower middle panel) and Helll 
(lower right panel) at times t/yr = [1 X 10 7 , 1 X 10* , 1 X 10 9 , 1 X 10 10 ] as indicated in the legend. We also show two equilibrium solutions from CLOUDY 
(with OTS approximation in green and full RT in blue). The thin black line is the C 2 -RAY result at t = 10 10 yr when using the U-OTS approximation. 
Note the larger difference between the two curves from C 2 -RAY after 10 10 yrs compared to the difference of the two curves from CLOUDY and the good 
agreement between C 2 -RAY using OTS with CLOUDY using OTS. 



without (black) and with (red) helium included. It can be seen that 
the power spectra are almost identical. We also show the relative 
difference, defined as loaW i^^^ H ^) '°» 10 ( A ij( H + fle )) global 

loglQ(&lJH + He))) ° 

ionization fraction of (x) ~ 0.1. The relative error is everywhere 
below 1%. This is similar for the other global ionization fractions. 
Therefore we conclude that the simplification of only using hydro- 
gen for reionization simulations with only stellar-type sources is 
legitimate. 



4.3.2 TEST 4: Multiple sources in a small cosmological density 
field 

This test was fully described in 1+06 and has subsequently 
been used in many pape r s on radiative transfer methods 
(e.g. ICantalupo & Porciaml 1201 ll : IPetkova & Springell 1201 ll ; 
Pawlik & Schavd boTUT Here we present results for this test 
including a cosmological helium abundance. We want to use this to 
illustrate the effect the presence of helium has on this test problem, 
as well as to compare to the original C 2 -RAY results from 1+06 
which used the simpler heating rate calculation, see Sect. 13.41 
Here we only summarize the most important aspects of the test 
setup and refer the reader to 1+06 for details: The density field is a 
snapshot at redshift z ~ 8.85. The simulation box has a side length 
of 0.5 /h comoving Mpc and the radiative transfer grid consists of 
128 3 uniform cells. The box boundaries are transmissive. The 16 
most massive halos in the box constitute the 16 sources. They have 
a constant photon output during the course of the simulation, the 



joint ionizing photon rate of all 16 sources is 3.29 x 10 ionizing 
photons per second. All sources have a black body spectrum 
with an effective temperature T c s = 100 000 K. The initial gas 
temperature in all cells is Tj n i = 100 K. All conditions are as in 
test4 of 1+06 except that our simulation has a helium abundance of 
nuc/n = 0.074. 

In Fig. [8] we show slices through the center of the simulation 
volume, just as shown in 1+06. We show the HI, Hel and Hell frac- 
tions (in a logarithmic colour scale) and the temperature (in a linear 
colour scale) together with the original C 2 -RAY results presented 
in 1+06 (figures 31 and 32) after 0.05 Myr of evolution. Before 
describing the visible differences in the results, we need to point 
out two additional (apart from the helium) important differences 
which mostly affect the heating: As described briefly in Sect. |3.4| 
the multi-frequency implementation of the heating forces us to use 
a timestep close to the ionization time scale. The original C 2 -RAY 
implementation used in this test did not have this restriction since it 
used the constant heating per photo-ionization approach. The C 2 - 
RAY results presented in 1+06 (right hand panels of Fig. [8]l used a 
timestep At = 0.001 Myr. We now use a timestep At = 0.00025 
Myr which was chosen on the basis of convergence studies. 

As can be seen in Fig. [8] the effect of using the multi- 
frequency heating instead of the constant energy- per-ionization 
heating is a lower temperature inside the HII region. Partly, this 
is also due to helium since the higher ionization energy of helium 
results in slightly less energetic photons. However, it can be seen 
that high density filaments close to the sources are in fact warmer 
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Figure 6. Results from TEST 2. Fractions of neutral/ionized hydrogen, neutral/single and double ionized helium and temperature as a function of distance 
to source for a source with /3 = 2 (red) and = 1 (blue) with temperature evolution and secondary ionizations for C 2 -RAY after t/yr = 10 5 , 10 7 and 10 9 
(decreasing line thickness). For comparison, we also include the equilibrium solutions of CLOUDY (dashed lines). 
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Figure 7. Results from TEST 3. Power spectrum of the ionized fraction x H + of simulation 53Mpc_g8.7_130S without (black lines) and with (red lines) 
helium included at four different global (mass averaged) hydrogen ionization fractions (x) as indicated in the legend. The differences are small: As can be 



seen in the top panel, for (x) ~ 0.1, the relative difference 
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are below 1 %. 
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than the less denser regions between them, while it was vice versa 
in the original implementation. This can be explained as follows: 
While in the low density region close to sources mainly low en- 
ergy photons with a higher ionization cross-section are absorbed, 
in the dense filaments, higher energy photons are absorbed which 
deposit more energy in the gas. However, dense knots in the fil- 
aments which do not host sources, still act as shields against the 
ionizing and heating radiation, resulting in embedded cold neutral 
regions. Another obvious difference is the less extended heating 
front outside the H II region. This is solely due to the inclusion 
of helium and its contribution to the optical depth, as tests without 
helium have shown. 

For the hydrogen ionization fraction field, we note that al- 
though the temperature inside the HII region is lower than in the 
original simulation, the ionization fractions are similar. This is due 
to the recombination photons from helium. In general, the differ- 
ences in hydrogen ionization fraction are very small everywhere. 
Due to the rather high effective temperature, T e g = 100 000 K, a 
considerable amount of photons capable of doubly ionizing helium 
is produced. The resulting He III regions can be seen as holes in the 
He II fraction in the lower left panel of Fig. [8] 

The average temperatures inside the H II regions are now (with 
the multi-frequency treatmen t of the heating) mo re similar to the 
temperatures from CRASH dMaselli et ai1l2003l) in 1+06 and re- 
sults of TRAPHIC presented in IPawlik & Schavd ( 1201 lh . How- 
ever, if those results also show higher temperatures for dense fila- 
ments inside the ionization front is not clear. 

Temperatures outside the ionization fronts cannot be com- 
pared since the results shown here are obtained with helium. The 
inclusion of helium effectively prevents a preheating ahead of the 
ionization front. 



5 CONCLUSIONS 

We presented an extension of the radiati ve transfer and photo- 
ionization code C 2 RAY, first introduced in Melle ma et all {2006b ) 
as a method for calculating non-equilibrium hydrogen photo- 
ionization. The new version treats the transfer of ionizing photons 
in step with time-dependent photo-ionization of both hydrogen and 
helium using the full OTS approximation, multi-frequency ioniza- 
tion and heating, as well as secondary ionizations. We described 
in detail the new elements to C 2 -RAY, such as the linearized so- 
lution of the set of hydrogen and helium rate equations using the 
coupled OTS treatment, and the calculation of the multi-frequency 
ionization and heating rates for both hydrogen and helium. 

We validated our implementation of the various ionization and 
heating processes, including the OTS approximation and secondary 
ionizations, by comparing to results from the photo-ionization equi- 
librium code CLOUDY. We validated our time-dependent solu- 
tions through convergence studies, which also provide us with 
timestep constraints. We confirmed that the new version of C 2 -RAY 
retains the property of being able to calculate ionization fractions 
with an accuracy of a few percent even for timesteps as large as 
0.1 times the recombination time and hard spectra. However, to ob- 
tain accurate temperatures, we found the timestep constraints to be 
more stringent. In the worst case timesteps of about 10% of the 
ionization time are needed to correctly determine the temperature. 
This condition can be relaxed in case the cells are optically thin, 
or when one is interested in temperatures beyond the recombina- 
tion time, or when the radiative cooling time is of the order of the 
timestep. 



A comparison between results obtained with the OTS approx- 
imation and with the often used, simpler U-OTS approximation 
(perhaps better known as the a B approximation), shows that there 
are significant differences in the results and that these differences 
are larger than those between the OTS approximation and the case 
of full radiative transfer of recombination photons. This implies 
that it is more accurate to use the full OTS approximation. We also 
verified that secondary ionizations have significant effects for hard 
enough spectra and should therefore be included. 

We used the new version of C 2 -RAY to re-run a cosmological 
reionization radiative transfer simulation with stellar-type sources 
that previously had been simulated with only hydrogen. By com- 
paring the power spectra of the HII fractions, we concluded that 
the morphologies of HII regions are not different between these 
two simulations. This confirms that it is correct to assume that for 
stellar-type sources, the helium ionization follows that of hydrogen. 

Finally, we presented the first version with helium for "Test 
4" from the Cosmological Radiative Transfer Comparison Project 
(1+06). We find that the inclusion of helium significantly affects the 
temperature distribution as the heating front is found to be steeper 
than in the hydrogen only case. This is caused by the larger cross- 
section of neutral and singly ionized helium at higher photon ener- 
gies. Helium is therefore stopping hard photons which would oth- 
erwise preheat the material far ahead of the the ionization front. 
This result confirms the original motivation of including helium in 
C 2 -RAY, it is an important absorber of EUV and SX photons and 
should therefore be taken into account when studying hard ionizing 
spectra. 

We are planning to use the new version of C 2 -RAY to explore 
the effects of quasar-like sources on reionization, both on the morh- 
phologies of HII regions and the heating of the IGM. However, the 
new version can be used to study any kind of photo-ionization prob- 
lem in which substantial amounts of Helll are formed, for example 
the inner parts of galactic HII regions around O and B stars. 
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APPENDIX A: DETAILED SOLUTION TO THE 
UNCOUPLED RATE EQUATIONS 

For completeness, we present here the solution to the helium rate 
equations in the uncoupled case (which we introduced as the U- 
OTS approximation in Sect. 13.1. It , i.e. the expressions for the 
eigenvalues, eigenvectors, particular solution vectors and coeffi- 
cients that can be used in the general solution to Eq. l(3j which is 
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Figure 8. Resuls from TEST 4. Slices through the center of the simulation box of TEST4 1+06 after 0.5 Myr. Upper panels from left to right: Hel, HI, HI 
from the original C 2 -RAY (for details see text); lower panels from left to right: Hell, temperature, temperature from the original C 2 -Ray. The color coding 
for H and He is logarithmic as indicated by the upper color scale, the temperature is in linear scale as indicated by the lower color scale on the right hand side. 



given by Eq. l[4}. For hydrogen, the values were given in Eq. ©. 
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This solution was also given in lAltav et alj J2008I) . 



I Schmidt ^Vb igt & Koeppen jl987t) who introduce the same 
general problem, Eq. l[3}, split up the matrix in two 2-state sys- 
tems and introduced a coupling between the stages which results in 
a single time dependence for the two stages of helium . This solu- 
tion w a s used for exam p le by Frank & Mellemal 1 1994 ), Raga et al] 
dl997l). IMellema et alj ( 11998) . iMellema & Lundqvistl d2002l) and 
IShapiro et alj <20041) . Here, we include explicitly the two different 
time dependences (the two different exponential constants), namely 
the two eigenvalues \Y e and \^ e . Their contribution is weighted 
differently for the two species of helium, see Eq. {4}. This is the 
mathematically correct solution and, as an extreme example sug- 
gests, also makes sense physically: Assume for the ionization rate 
ThoU — > and initially, thcIii > 0, say XHeiii = 0.5. In this 
case, the time evolution of iHciii should not (directly) depend on 
the ionization rate ThcI, so the time evolution of shcII and xnein 
should not have the same exponential factor. 



The averages of the ionization fractions over one timestep At 
can be calculated for both hydrogen and for helium as: 

, 4 fx(t)dt A , „, r^At_^ 

(as) = f ]_ jM = Ci*i 21, + x p with 2U = 



/ tdt 



At Ai 
(All) 



To avoid floating point precision problems is set to 1 explicitly 
when XiAt < 10" 8 . 
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APPENDIX B: DETAILED SOLUTION TO THE 
COUPLED RATE EQUATIONS 

In the case of the OTS approximation, the rate equations for hy- 
drogen and helium are coupled by the recombination photons from 
the different ionization stages of helium. The solution to Eq. {3} 
with the elements defined as in Eq. j 1 2t and Eq. d 1 3 1 > is given by 
Eq. © with the following values for the eigenvalues, eigenvectors, 
particular solution vector and coefficients a : 
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Here, the coefficients S, K, R and T are defined as 

S = 
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The averages can be calculated with Eq. dAl 1 b where Ai, a, x p and 
Xi as introduced in this Appendix. 



APPENDIX C: IONIZATION CROSS SECTIONS 



As described in Section [331 within each sub-bin we assume the 
same frequency dependence of the ionization cross sections for 
all three species. In all sub-bins in bin 2, we use the power law 
indices fro m our fit to the ne utral helium ionization cross-section 
data from lVerner et alj q 19961) . In the sub-bins from bin 3, we use 
the power law indices from the fit to the ionized helium ionization 
cross-section. To fit the power-law, we use a linear-least square fit in 
the (log 10 v, log 10 a) space. The power-law indices are presented 
in fi sure [CTl for (1,26,20) subbins in bin (1,2,3). 
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Figure Cl. Power law indices of HI (b lue), Hel (red) and Hell (green) 
ionization cross-section from the fits to IVerner et alj (l996). The vertical 
lines indicate the threshold frequencies. Note the very different power law 
indices for HI and He I in frequency bin 2 and the very similar power law 
indices for all three species above 10 17 Hz. 



APPENDIX D: ON THE DIVISION OF PHOTONS 

Several different suggestions for the distribution of the ionizing 
photons between different species have been proposed in the lit- 
erature. We follo w the method pr oposed by A. Lidz (private com- 
munication) and lAltav et alj | |2008|) . who distribute the photons de- 
pending on the ratios of the optical depth. For two species labelled 
1 and 2: 



n 



(Dl) 



Tl + T'2 

where A is the fraction of io nizations going into ionizing species 1 
with n . iBolton et alj (2004) suggested the following recipe 



h 



(1 - cxp(— Ti))exp(-T 2 ) 



(1 - exp(-ri)) exp(-r 2 ) + (1 - cxp(-r 2 )) cxp(-n) 



and lMaselli et alj j2003h chose 

1 - exp(-n) 



h 



(1 - exp(-Ti)) + (1 - exp(-T 2 )) 



(D2) 



(D3) 



By considering the limits of low-optical depth limit and high 
optical depth limit it becomes clear that Eq. JDlb is the correct one: 
For the low optical depth limit, the intensity leaving a cell with op- 
tical depth t is I(t) = Jo (1— r), where Io is the ingoing intensity. 
The absorbed fraction, (Io—I(t))/Iq depends therefore linearly on 
r. For very high optical depths, the fractions according to Eq. (D21 
and Eq. dD3t become largely independent of the fraction of actual 
optical depth of the two species: Independent of the exact ratio of 
the optical depth, Eq. jD3\ yields values close to 0.5 for both frac- 
tions while the fractions are close to 1 and close to for Eq. dD2l l. 
This cannot be correct as the fractions should still depend on the 
relative values of r\ and r 2 . This is illustrated in Fig. lDll where we 
plot the fraction A as a function of optical depth of both species 
for all three recipes. 

We would like to point out that the latest version of CRASH 
now also uses Eq. l IDU to distribute the photo-ionizations over the 
different species (Maselli, private communication). 
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Figure Dl. For a two species medium, we show the fraction of ionizing 
photons absorbed by species 1 as a function of the optical depth of both 
species using the three different recipes, L (Eq.|DT] left panel), C (Eq.|D3l 
middle panel) and B (Eq. lD2l right panel). The fraction is color coded as 
indicated by the color bar, red corresponding to /i =1 and blue corre- 
sponding to fi =0. 

APPENDIX E: VARYING THE NUMBER OF SUB-BINS IN 
THE FREQUENCY BINS 2 AND 3 

As mentioned in Section [3~3l the code can use different numbers of 
sub-bins in frequency bins 2 and 3. As described in Appendix [Cl 
the code uses a single power law fit for all species in each sub-bin. 

This has the following disadvantage: In frequency bin 2, for 
example, a single power law fit to the ionization cross-sections of 
each species (separately), He I and H I would be a sufficiently good 
fit and yield power law indices: Shi = 2.91 and sh c i = 1.88. 
However, since these are very different numbers, using 1.88 as a 
power law index over the entire bin 2 for both species, introduces a 
substantial error in hydrogen optical depth. Therefore, we introduce 
sub-bins. 

We tested different numbers of sub-bins for bins 2 and 3. For 
this test, we use the parameters from TEST2 with f3 = in or- 
der to have a substantial fraction of ionizing photons in bin 3. We 
show our results in Fig. IE 1 1 where the four panels on the left hand 
side show the relative differences of the ionization fractions (H II, 
He I, He II and He III) using 1,2,3,6 and 10 sub-bins in bin 2 as 
compared to using 26 sub-bins. The four panels on the right hand 
side show the relative difference of the ionization fractions using 
1,4,9,1 1 and 16 sub-bins in bin 3 as compared to using 20 sub-bins. 
For bin 2, it can be seen that the relative differences using 10 sub- 
bins compared to 26 is at most at 4 %. This maximum error is at 
the ionization front position. Here, the error introduced by using 
the OTS approximation as opposed to including diffuse photons is 
most probably higher than that. Therefore, we conclude that using 
more than 10 sub-bins in bin 2 might not sufficiently improve the 
results. For bin 3, it can be seen that the relative differences using 
1 1 sub-bins compared to 20 is at most 2 % which is reached far 
outside the front in the H II and He III fractions. Apart from those 
locations, the error is well below the percentage level. We there- 
fore conclude that for most applications, 1 1 sub-bins in bin 2 are 
sufficient. 



APPENDIX F: DEPENDENCE ON TIMESTEP 

As was pointed out in M+06, the approximation of using time- 
averaged values for the ionization states in the calculation of the 
ionization rates to avoid the need of small timesteps is strictly only 
valid in the case of negligible contribution from collisional ion- 
izations and recombinations. For those processes it does matter at 
which time during the timestep, they occur. 

Given the added complexity due to the inclusion of helium, 



the coupled OTS approximation and the multi-frequency photo- 
heating, we present in this Appendix convergence tests for the 
timestep. We first consider the convergence of the ionization frac- 
tions at constant temperature and then the convergence of the tem- 
perature evolution. 

We test the effect on the ionization fractions of varying the 
timestep 5 orders of magnitude for a source with a power-law spec- 
trum with power law index /3 = 1. Fig. lFll shows the relative error 
for a test with the same parameters as in TEST2 at t = 10 s yr 
(this corresponds roughly to the recombination time scale) using 
At = 10 5 , 10 6 , 10 7 and 10 s yr, compared to At = 10 3 . The lat- 
ter corresponds roughly to the ionization time for the first cell. As 
can be seen, even for the large timestep At = 10 7 , the maximum 
error is 3% and for At = 10 6 yr, the error is everywhere well be- 
low the percent level. We therefore conclude that At < 0.1t ICC is a 
sufficient timestep criterion for accurate ionization calculations. 

Next, we test the temperature evolution dependence on the 
timestep. Here, we show the results of two one dimensional sim- 
ulations, one with optically thin cells (r ~ 0.6) and one with mod- 
erately optically thick cells (r ~ 60). Both have recombination 
timescales t ICC ~ 10 5 yr. The ionization time scale for the first cell 
for the simulation with optical thin cells is ti on ~ 5 x 10~ 3 yr 
and the the moderately optically thick case has tion ~ 0.5 yr. The 
parameters for these tests are nu = 0.926 cm -3 , uhc = 0.074 
cm -3 ,7V 7 = 10 46 (10 48 ), Ar = 10 12 (10 14 ) km, T eS = 100000 
K, Tlni = 100 K for the optically thin (moderately thick) case. We 
use (ri2,n3)=(26,20) frequency sub-bins. 

In Fig. lF2l we present the results for these t wo tests. The form 
was in spired by that of testO in 1+06 and testl in Paw lik & Schavd 
fcOllh . For four cells of our computational grid we show the 
temperature evolution for six different choices of timestep: At/yr 
= 10~ 5 , 10" 3 , 10~\ 10 1 , 10 3 and 10 5 . As we evolve each case 
only for 10 J timesteps, the results of each simulation only overlaps 
with two others. When two curves overlap only the one with the 
longest timestep is shown. The lower rows show the relative error 
between two subsequent choices of timestep. 

We see from this figure that while timesteps larger than ti on 
still yield reasonable good results for the optically thin case, this is 
not true for the optically thick case. For the optically thin case, the 
error increases with distance to the source. This is due to an under- 
estimation of the preheating of cells further away from the source. 
For those the heating rate as function of time is already a substan- 
tial fraction of its maximum before the front reaches the cell. This 
dependence on the optical depth between a cell and the source pre- 
cludes "local" fixes based on the cell properties and evolution. 

For the optically thick case, timesteps of the order of ion- 
ization timescale result in errors of the order of 8% (for the first 
cell) compared to timescales several (2 and 4) orders of magnitude 
smaller than the ionization timescale. Although this error is larger, 
it actually decreases with distance, mostly because less preheating 
occurs at larger distances. 

If one is only interested in time-scales larger than the recom- 
bination time-scale, the errors are small in both the optically thick 
and optically thin case. This is most probably because the initial 
spike in the heating rate becomes a less important contribution to 
the total photon energy input in each cell. 
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Figure El. Relative difference of HII, Hel, Hell and Helll fractions. The input parameters are as in TEST1B. In the four panels to the left, the number of 
sub-bins in frequency-bin 2 is varied according to the legend. In the four panels on the right hand side, the number of sub-bins in frequency-bin 3 is varied 
according to the legend. 




Figure Fl. Relative difference of HII, Hel, Hell and Helll fractions at t = 10 8 yr. The timestep At was varied according to the legend. The comparison is 
made against a timestep that correspond roughly to the ionization time of the first cell, At = 10 3 yr. The input parameters are as in TEST2 with = 1 but 
without temperature evolution. The insets shows the ±1% regions. 



APPENDIX G: RECOMBINATION- COOLING AND 
COLLISIONAL IONIZATION RATES 



RECOMBINATION For the hydrogen recombination rates, a^n 
and a ml and for the helium recombinat ion rates, Q^hpttt an d 
a Hciii' we are using the fitting formula from Hui & Gnedinlll997h . 
hencefort h HG9 7. These are also good fits to the data from 
iHummed i 1994) (accurate to 1.5 % in the temperature range 



prensented there, 10 K to 10 K) even if t hey were fitted to the 
slightly older data from lFerland et alj jl992h . 

The recombination coefficients OhoII an d "Hen above T = 
7 x 10 4 K are dominated by dielectronic recombination. We in- 
clude it above T = 1.5 x 10 4 K acc ording to the fitting formula 
from lAldrovandi & PequignoJ d 19731) as an extra contribution to 
the recombination rates. Below T = 9 x 10 3 K, the fitting formula 
for Qhii an d ami from HG97 provide a good fit (accurate to 6 % 
from T = 10 K to 9 x 10 3 K) to the data for a£ eII and afeii 
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Figure F2. Upper panel: Temporal evolution of the temperature for 4 optically thin (t ~ 0.6) cells (panels left to right correspond to cells 1, 4, 7 and 10). 
Lower panel: the same for moderately optical thick cells (t ~ 60). We show results for 6 different timesteps according to the legend. All times are in years. 
With each timestep At we evolved the box for a time corresponding to 10 5 X At. We also show the relative error between the results using two consecutive 
timestep sizes, according to the legend: relative error = T ( At ~ H*gj_Z^i^= 10 y r ) Indicated on the abscissa (and by the dotted vertical grid lines) are the 
recombination time scale and the ionization time scale (of the first cell). Note that the relative error for the optically thin case is larger for cells further away of 
the source while the error decreases with distance to the source for the optically thick case. 
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from lHummer & Storevl Jl998h . Above T = 9 x 10 3 K, we are 
using the fits for QhoII an d q h c ii from HG97 which provide fits 
accurate to 6 % up to T = 2.5 x 10 4 K, the highest temperature 
Hum mer & Storevl i 19981) provides data for. 

For th e recombination to n = 2 , o^eiii, we fit me 
data from lOsterbrock & Ferland (2006) as a HcIII (T) = 

3.4 x io- 13 (r/io 4 A')- - 6 . 

COLLISIONAL IONIZATION For the collisional ionization 
coefficients Chi, ChcI and ChcII we use the fitting formulas from 
HG97 in the temperature ra nge 2 x 10 5 K < T < 10 7 K where 
Chi (Chgi) fit the data from ljanev et alj i 19871) to an accura cy of 
12% (5 %). Below 2 x 10 s we use the fitting formula from lCoxl 
( 1970) which is simpler and slightly more accurate than HG97 for 
lower temperatures. 

COOLING For cooling coefficients for HII and Helll (free-free 
+ reco mbination cooling) we interpolate the data from iHumme i 
(1994). For the HI cooling coefficient we inclu de collisional ex - 
citation cooling and use the line strength from lAggarwall Jl983l) . 
For the Hell cooling coefficient we include collisional excita- 
tion, collisional ionization and dielectronic recombination cool- 
ing from HG97 and B re combination and free-free cooling from 
Hummer & Storevl J 1998h . For Hel we only include collisional ion- 
ization (from HG97) since its col lisional exci tation is negligible 
and the rate proportional to n e (c.f. lBlacklll98lh . 

In case of cosmologi cal simulations, we include Compton 
cooling according to Shapi ro~& Kand Jl987h as well as cooling due 
to cosmological expansion. 
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